#include "tmatrix.h"
#include <stdlib.h>


Tmatrix::Tmatrix() {
    a11.rl=0;a12.rl=0;a21.rl=0;a22.rl=0;
    a11.im=0;a12.im=0;a21.im=0;a22.im=0;
}

Tmatrix::Tmatrix(Complex k1,Complex k2,double m1,double m2,double d) {
    Complex z, i(0.0,1.0), tmp;
    Complex coeff;
    z	=  (m2/m1)*(k1/k2);
    coeff = 0.5*cplx_sqrt((m1/m2)*(k2/k1));
    a11	= coeff*(1+z)*cplx_exp(i*(k1-k2)*d);
    a12	= coeff*(1-z)*cplx_exp(i*(k1+k2)*(-d));
    a21	= coeff*(1-z)*cplx_exp(i*(k1+k2)*d);
    a22	= coeff*(1+z)*cplx_exp(i*(k1-k2)*(-d));
    if(	a11.rl==0 && a11.im==0 &&a12.rl==0 && a12.im==0 && 
	a22.rl==0 && a22.im==0 && a21.rl==0 && a21.im==0) {
      printf("TMATRIX: Zero matrix\n");
      k1.show();
      k2.show();
      exit(1);
    }
    tmp = a11*a22-a12*a21;
    if(tmp.rl==0 && tmp.im==0) {
      printf("TMATRIX: Matrix with zero determinant\n");
      show();
      k1.show();
      k2.show();
      printf("m1=%lf, m2=%lf, d=%le\n",m1,m2,d);
      exit(1);
    }
}

Tmatrix Tmatrix::id() {
    Complex unit(1.0,0.0),null(0.0,0.0);
    a11	= unit;
    a12	= null;
    a21	= null;
    a22	= unit;
    return *this;
}
Tmatrix Tmatrix::operator=(const Tmatrix & op)	{		
    a11	= op.a11;
    a12	= op.a12;
    a21	= op.a21;
    a22	= op.a22;
    return *this;
}
Tmatrix  Tmatrix::operator*(const Tmatrix & op)	{
    Tmatrix	temp;
    temp.a11	= a11*op.a11+a12*op.a21;
    temp.a12	= a11*op.a12+a12*op.a22;
    temp.a21	= a21*op.a11+a22*op.a21;
    temp.a22	= a21*op.a12+a22*op.a22;
    return	temp;
}
Tmatrix Tmatrix::is(Complex k1, Complex k2, double m1, double m2, double d)    {

    Complex z,i(0.0,1.0),tmp;
    Complex  coeff;

    if(k2.rl==0.0 && k2.im==0.0 ) {
      printf("TMATRIX.IS(): ERROR ... k2==0\n");
      printf(" k1=");
      k1.show();
      printf(" and d= %le\n",d);
      exit(1);
    }
    if(k1.rl==0.0 && k1.im==0.0 ) {
      printf("TMATRIX.IS(): ERROR ... k1==0\n");
      printf(" k2=");
      k2.show();
      printf(" and d= %le\n",d);
      exit(1);
    }
    z   =  (m2/m1)*(k1/k2);
    coeff = 0.5*cplx_sqrt((m1/m2)*(k2/k1));
    a11 = coeff*(1+z)*cplx_exp(i*(k1-k2)*d);
    a12 = coeff*(1-z)*cplx_exp(i*(k1+k2)*(-d));
    a21 = coeff*(1-z)*cplx_exp(i*(k1+k2)*d);
    a22 = coeff*(1+z)*cplx_exp(i*(k1-k2)*(-d));

    if( a11.rl==0 && a11.im==0 
        && a12.rl==0 && a12.im==0 
        && a22.rl==0 && a22.im==0 
        && a21.rl==0 && a21.im==0) {
      printf("TMATRIX.IS(): Zero matrix\n");
      k1.show();
      k2.show();
      exit(1);
    }
    tmp = a11*a22-a12*a21;
    if(tmp.rl==0 && tmp.im==0) {
      printf("TMATRIX.IS(): Matrix with zero determinant\n");
      show();
      k1.show();
      k2.show();
      exit(1);
    }
    return *this;
}

void Tmatrix::show(void) {	
  printf("T-MATRIX:\n");
  printf("---------------------------------------------------------------------\n");
  printf("| %le + i %le   |   %le + i %le |\n",a11.rl,a11.im,a12.rl,a12.im);
  printf("---------------------------------------------------------------------\n");
  printf("| %le + i %le   |   %le + i %le |\n",a21.rl,a21.im,a22.rl,a22.im);
  printf("---------------------------------------------------------------------\n");
  printf("\n");
}

void Tmatrix::show(FILE * fp) {	
  fprintf(fp,"T-MATRIX:\n");
  fprintf(fp,"---------------------------------------------------------------------\n");
  fprintf(fp,"| %le + i %le   |   %le + i %le |\n",a11.rl,a11.im,a12.rl,a12.im);
  fprintf(fp,"---------------------------------------------------------------------\n");
  fprintf(fp,"| %le + i %le   |   %le + i %le |\n",a21.rl,a21.im,a22.rl,a22.im);
  fprintf(fp,"---------------------------------------------------------------------\n");
  fprintf(fp,"\n");
}


/* ---------------------------------------------------------------- */
/* ---------------------------------------------------------------- */
Complex	product(int num, Tmatrix T[])	{

  int j;
  Tmatrix temp;
  Complex tmp,tmp2;

  if(num>1000)	{
    printf("PRODUCT: trouble. Too many T-matrices! Increase the limit in file tmatrix.cpp and recompile\n");
    exit(1);
  }
  temp = T[0];
  for(j=1;j<num;j++)	{
    temp = T[j]*temp;
    tmp2 = temp.a11*temp.a22-temp.a12*temp.a21;   
  }
  if(temp.a22.rl==0.0 && temp.a22.im==0.0) {
    printf("PRODUCT: Tmatrix with a22 == 0 \n");
    for(j=0;j<num;j++) T[j].show();
    temp.show();
    exit(1);
  }
  tmp = (temp.a11*temp.a22-temp.a12*temp.a21)/temp.a22;
  return tmp;

} 









